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ABSTRACT: We consider the detection and initial guess problems for the LISA gravitational 
wave detector. The detection problem is the problem of how to determine if there is a signal 
present in instrumental data and how to identify it. Because of the Doppler and plane- 
precession spreading of the spectral power of the LISA signal, the usual power spectrum 
approach to detection will have difficulty identifying sources. A better method must be 
found. The initial guess problem involves how to generate o priori values for the parameters 
of a parameter-estimation problem that are close enough to the final values for a linear least- 
squares estimator to converge to the correct result. A useful approach to simultaneously 
solving the detection and initial guess problems for LISA is to divide the sky into many pixels 
and to demodulate the Doppler spreading for each set of pixel coordinates. The demodulated 
power spectra may then be searched for spectral features. We demonstrate that the procedure 
works well as a first step in the search for gravitational waves from monochromatic binaries. 



I. The detection problem 

The signal coming down from a spaccbornc gravitational wave detector will contain 
signals from many close compact binaries, but, to look at the signal, one would see nothing 
but noise. This is because the instrument noise, integrated over the bandwidth of the detector, 
far exceeds the strength of any individual binary signal. This situation is common in data 
analysis, and the usual procedure for detecting weak coherent signals in the presence of 
large incoherent noise is to generate the power spectrum of the signal. Since the binary 
star signal is monochromatic, its power will be concentrated in a single frequency bin of the 
power spectrum, while the noise is spread over all frequencies. The signal will then be seen 
as a single spike rising above the neighboring noise spectrum. However, the usefulness of 
this power spectrum method for signal detection is limited for the spaceborne gravitational 
wave detectors, due to the motion of the detector around the sun and to the change in the 
orientation of the detector. As a result of these motions, the gravitational wave signal from 
a particular source will be strongly phase modulated, spreading the total power from the 
source over many frequency bins. It is therefore quite possible that none of the frequency 
components of the signal will be seen above the noise. 

To illustrate, we show in Fig. 1 the power spectrum of the gravitational wave signal from 
a binary star, where the signal frequency is 0.01 Hz and the amplitude is 3.5 x 10~^^. The 
signal has been sampled at 10 s intervals for one year. This is the power spectrum that would 
be seen if the detector were at rest in the solar system. The power spectrum as seen in a 
moving detector is shown in Fig. 2. To generate this figure, the detector has been assumed to 
be the LISA detector [1], in which three free-flying spacecraft move on heliocentric trajectories 
in such a way as to remain at roughly constant distances from each other and to remain in 
a plane inclined by 60° to the plane of the ecliptic, the plane precessing about the ecliptic 
pole once per year. The spreading of the power over about 100 frequency bins can be clearly 
seen in the figure. If a noise spectrum, characteristic of the proposed LISA instrument noise, 
is added to the signal, the power spectrum in Fig. 3 results. As may be seen, the spreading 
of the binary signal power over many frequency bins has produced a power spectrum that 
cannot now be seen above the noise. 

One of the goals of this paper is to show how this "detection problem" may be solved. 

II. The initial guess problem 

The ultimate goal of gravitational wave detectors is to use the gravitational wave sig- 
nals to determine the astrophysical properties of the gravitational wave source. Papers by 
Cutler and Vecchio [2,3] and by Moore and Hellings [4] have addressed the problem of what 
information is available in the signal and how well the parameters of the source may be de- 
termined from a signal with a particular signal-to-noise ratio. The data analysis technique 
to be used to extract the information from the signal is called linear least-squares parameter 
estimation. This technique accounts for the possible correlation of the parameters with each 
other and gives the best-fit values of the parameters and their formal standard deviation, 
with the assumption of a background of Gaussian noise. However, the "linear" in "linear 
least-squares" reminds us that this procedure will only work in a totally linear problem or in 
the limit that the problem may be linearized by writing it in terms of small quantities, the 
differences between the final best-fit values for the parameters and some a priori values. 

The signal from a monochromatic binary is characterized by seven independent param- 



eters. These are: 



• the frequency / of the gravitational wave signal (a parameter that depends on the 
masses of the two members of the binary and on the radius of their mutual orbit) , 

• the amplitude h (which depends of the binary masses, the orbital radius, and on the 

distance from the binary system to the Earth), 

• the two-parameter location of the binary in the sky (ecliptic coordinates 9 and (j)), 

• the two-parameter orientation of the orbit plane of the binary (taken to be the incli- 
nation i of the orbit plane normal to the line-of-sight to the binary and the angle if) 
by which the major axis of the binary's apparent ellipse in the sky is rotated relative 
to the ecliptic plane, around the line-of-sight), and 

• the phase 0o of the binary orbit at some initial time to = 0. 

The received gravitational wave signal is linear in only one of these seven parameters, the 
amplitude h. For all the rest, initial guesses must be found that are sufficiently close to the 
correct values of the parameters to allow the linear least-squares analysis to converge to the 
proper place. Often, in parameter estimation processes like this, the initial guesses are known 
from some other aspect of the observations. For example, in a similar problem to the one we 
are facing, the times-of-arrival of pulsar pulses are used to determine the position of the pulsar 
in the sky and the parameters characterizing pulsar's rotational dynamics. The initial guess 
for the position, however, is not found from the time series itself, but from the direction in 
which the antenna is pointing when the pulsar data is collected. The spaceborne gravitational 
wave detector is an omnidirectional instrument, so there will be no such ancillary information 
on source direction; all information must be extracted from the time series itself. 

A critical question for LISA data analysis, therefore, is how does one take the time series 
and generate initial guesses for the six non-linear parameters. One method would be to cross- 
correlate the time series with a series of templates that represent all possible values of all 
possible parameters and look for the templates that have a strong correlation with the signal, 
the parameters of the template becoming the initial guess. The spacing between templates 
is set by the need to get an initial guess that is within the linear regime of the final answer. 
Ultimately, this is the only method that can be employed in the absence of data from outside 
the time series. However, let us consider the nature of the six non-linear parameters to see if 
a more step-by-step method may be found. 

The frequency / of a signal buried in noise is ordinarily found by spectral analysis. The 
reason why that avenue is not available here is evident by inspection of Figs. 2 and 3; the 
spectral power is spread out over many frequency bins. The contribution to this spreading 
by the motion of the detector around the sun is given by the usual Doppler formula 



fr = f 



1 sin 9 sin(r2i — 0) 



(1) 



where is the frequency observed in the detector, / is the frequency that would be observed 
at the solar system bary center, Q is the mean motion of the detector around the sun, and 
R is the radius of the detector orbit. The magnitude of the Doppler spreading depends on 
the gravitational wave frequency / and on 9. The phase of the Doppler depends on 0. For 
gravitational waves at frequencies higher than / = c/{2R){k, 10~^ Hz for LISA) this Doppler 
shift dominates over the spreading produced by the precession of the detector plane and is 



the major impediment to determination of the frequency by spectral analysis. Now, if the 
frequency and position of a source were to be known, the template matching could be a much 
simpler procedure, with only three parameters left to determine -1,1^, and ^o- What is 
more, these three parameters have little effect on the signal in the detector and are much less 
well-determined in the final least-squares analysis* than are 9 and (p. Therefore the resolution 
needed to find initial values for these parameters that lie within the linear regime, where next- 
higher-order terms are not greater than the final uncertainty, is very coarse, and the number 
of templates needed to span the parameter space will be small. A method to independently 
find the position and frequency would thus greatly simplify the template matching procedure. 

The second goal of this paper is to show how to generate such a step-by-stcp approach 
to the "initial guess" problem by approximating the frequency and location of the source 
independently of the other parameters of the system. 



III. Demodulating the detector output 

Both problems we have discussed, the detection and initial guess problems, may be 
solved by first undoing the position-dependent Doppler modulation the signal undergoes. The 
Doppler modulation depends only on 9 and (p and is given simply by Eq. 1. The modulation 
produced by the detector-plane precession is much more complicated (see Moore and Hellings 
[4]) and depends on all four angular parameters, 9, (f), i, and ip, so demodulating the Doppler 
only is both a simple task to undertake and one that will undo the greatest part of the 
spectral spreading at frequencies where the spreading is most pronounced. 

Given a position on the sky, and given a time series from the detector, a new time 
series is generated which is the time series that would be received from that direction if the 
detector were at rest at the solar system barycenter. The plane of this fictitious detector will 
continue to precess, however, at the yearly period, just like the true LISA detector. If the 
time series at the detector is y(t„), we seek a time series z{tn) that represents the signal that 
the barycentric detector would have received. The algorithm for generating z{tn) is 

Z{tn)=y{tn-T), (2) 

where r is the time of flight from when the wavefront passed the detector to when it arrives 
at the barycenter, as seen in Fig. 4. The value of r, given comes from the solution of the 
transcendental equation 

r=^:%^, (3) 

where h is the unit vector toward the source and f{t) is the vector from the barycenter to the 
detector at time t. At time tn, the value of r is found by solving Eq. 3, and the barycentric 
time series is generated using Eq. 2. 

This procedure, of course, assumes that one knows the direction n to the source. When 
the direction is known, and the procedure applied, a time series is generated whose Fourier 

* For a typical case in the binary star solutions from Moore and Hellings [4], the uncer- 
tainties in radians in the angular parameters of a binary star with frequency 0.01 Hz and 
amplitude h = 2.2 x 10~^^ were ae = 3.68 x 10~^ and cr^ = 1.18 x 10~^ for the position 
parameters and (jj = 2.00 x 10^ and = 1.18 x 10^ (i.e., not determined) for the source 
orientation parameters. 



spectrum is shown in Fig. 5. As may be seen in comparison with Fig. 3, the spread spectrum 
that lay beneath the noise in the uncorrected time series is now seen strongly above the noise. 
The breadth of the spectral "line" in Fig. 5 is a result of the uncorrected modulation due to 
the precession of the detector plane. Indeed, when we produce a demodulated spectrum for 
the OMEGA detector [4], where the detector plane does not precess, the entire signal power 
is found in a single line, as seen in Fig. 6. 

Clearly, the remaining problem is that of determining the direction to the source. With 
no a priori knowledge, the only course of action is to consider all directions one at a time, to 
demodulate for that direction, and to examine the power spectrum of the resulting demodu- 
lated time scries. This has the effect of solving the detection problem and the initial guess 
problem simultaneously. The demodulated signals will collect the power into a few frequency 
bins where they will protrude above the noise spectral power, the frequency at which this 
occurs will be the frequency of the source, and the particular direction that produces this 
visible spectral feature will be the direction to the source. 

To illustrate this procedure we have simulated signals from three binary systems, as they 
would be observed in the LISA detector. We have combined these signals together, as they 
would appear in the detector, and then added random Gaussian noise at a level consistent 
with the LISA detector in this frequency band. The time series is sampled every 10 s and 
accumulated for a year (3.16 x 10^ data points). The sky is then divided into pixels of one 
square degree, for a total of 20,582 pixels. The pixel coordinates are used one at a time 
to demodulate the data according to Eqs. 2-3, and the power spectrum is formed for each 
case by fast Fourier transform of the demodulated time series. The power spectrum is then 
searched for features that are a chosen factor above the average power in the spectrum. Since 
the spectral lines are still spread over several adjoining bins by the detector plane precession, 
the detection procedure takes three-point running averages of the power spectrum to enhance 
such features as it searches. The result of this procedure is seen in Fig. 7. 

IV. Discussion 

The limits of the sky plot in Fig. 7 are from the ecliptic plane {6 = 90°) to the ecliptic 
pole = 0°) and over 360° in ecliptic longitude. The negative ecliptic lattitudes {9 > 90°) 
are the mirror image of the positive ones due to the symmetry of sin 9 in Eq. 1. The resolution 
of the plot is one square degree. The three sources, at {9,(/)} coordinates {60,120}, {78,20} 
and {12, 280} are seen clearly above the noise in the plots. However, the size and complicated 
structure of each region make it difficult to select a value for the coordinates of the source 
to the same accuracy as the resolution of the plot. In each of the three regions, a search for 
the brightest pixel for each source will give values that are within ±10° of the correct values, 
but this is not sufficient to be an initial guess for 9 and (j). This may be seen by considering 
the Doppler formula in Eq. 1. The change in the Doppler shift of the signal due to a change 
in 9, for example, is proportional to 

6(sm 9) = cos 9 69 — — sin 9 59^ + Higher Order Terms (4) 

The second term on the right represents the non-linear error introduced into the linear least- 
squares fit. For the source at 9 = 60°, this non-linear term would be 1.3 x 10~^ radians 
for 69 = 10°. This should be compared with the error in the first term due to the final 
uncertainty in the linear least-squares parameter fit, cos9ae = (0.5)(5 x 10~^) = 2.5 x 10~^ 



(from the least-squares fits of Moore and Hellings [4]). Thus, when a hnear least-squares fit 
tries to make a 10° change from initial guess to final value for 9, it will use only the first term 
in Eq. 4 and will make an error five times its calculated formal uncertainty. The situation for 
(f> is identical. Prom this, we conclude that one cannot go directly from the position estimates 
taken from the demodulated signal to a final least-squares fit. 

So what is to be done? At this point we need to remember that our goal in the de- 
modulation procedure was not to generate initial guesses directly, but to reduce the range of 
templates that would be needed for a template-matching procedure. Indeed, the complicated 
structure that keeps us from determining a single accurate position via the demodulation is 
a result of the precession of the detector plane that depends on 9, (j), i, and ip. For the non- 
precessing OMEGA detector, the demodulation alone is sufficient to get the initial position 
guess to the required accuracy, as may be seen by examination of Fig. 8. For LISA, the result 
of the demodulation has only been to find an exact value for the frequency and approximate 
values for 9 and 0. However, this will enable one to reduce the required template space from 
(10^ frequencies) X (2 x IC^ positions) to less than 100 position templates. A five-parameter 
template matching procedure (/ is now known and h is always linear) over a restricted set 
of values for 9 and 4> will be needed to produce a sufficient best guess to be supplied to the 
final least-squares estimator. 

Several other things may also be noted in Fig. 7. First, a single bright source will 
dominate the observed demodulated power over a large range of values of 9 and (p, possibly 
masking weaker nearby sources. Part of this is an artifact of the method we have used to 
generate Fig. 7, the picking out of the strongest spectral feature to represent each pixel. 
Perhaps the wisest use of the figure is to eliminate regions where there is no significant 
demodulated power as being a region free of sources that reach above the noise at all. One 
may then concentrate on the regions with significant power and examine spectra from many 
points in the region to look for secondary sources with different frequencies. A second thing 
to be noticed Fig. 7 is the significant leaking of power from sources near the ecliptic plane 
{9 = 90°) down towards the ecliptic plane. This is a result of the diminished position accuracy 
for the least squares fits for monochromatic binaries as one approaches the ecliptic (see Moore 
and Hellings [4], Fig. 3). Looking at it another way, this is a result of the fact that the phase 
of the Doppler modulation is the same for sources at the same ecliptic longitudes; only the 
amplitude of the modulation is different. 

Finally, let us consider some of the questions that remain to be addressed. First, the 
details of the final template matching need to be studied. Cornish and Larson [5] have 
generated an efficient code for template matching that works in the frequency domain and 
that is capable of fast correlation of signals over the entire parameter space. They have also 
generated a Doppler demodulation code[6], that similarly works in the frequency domain, 
whose efficiency should be compared with the time-domain code that we have used in this 
paper. Lastly, we should point out that among the test cases we have run are some in which 
there are two or more sources nearby each other in frequency and in position. When this 
occurs, there seems to be a pulling of the best demodulated pixels away from both source 
positions to some intermediate position. This has implications not only for the initial guess 
problem, but also for the problem of parameter estimation in the presence of non-Gaussian 
noise. These are issues we intend to address in future papers. 
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CAPTIONS: 



Figure 1. A portion of the power spectrum of the gravitational wave from a binary star, as seen 
by a stationary, non-preeessing deteetor at the baryeenter. The parameters characterizing 
the gravitational wave source are /i = 3.5 x 10~^^, / = 1.0 x 10~^ Hz, 9 = 60°, ^ = 120°, 
i = 10°, i/; = 84°, and = 0. 

Figure 2. The same gravitational wave as in Fig. 1, as seen by the LISA detector. The power 
is now spread over about one hundred frequency bins. 

Figure 3. The same gravitational wave in the LISA detector as in Fig. 2, but with noise added 
at a level characteristic of the LISA detector in this frequency band. The spread spectrum 

cannot now be detected above the noise. 

Figure 4. The geometry of the revision of the wavefront timing to the solar system baryeenter. 
In the first picture, the LISA detector receives the signal at time t = tn — r. In the second 
picture, a time r later, the same signal is received at the baryeenter, while LISA has moved 
to a new position. 

Figure 5. The signal from Fig. 3, after the revision to the solar system baryeenter has 
been performed for the particular direction of this source. The signal at 0.01 Hz is clearly 
visible above the noise. The remaining breadth of the spectral peak is due to the modulation 
produced by the precessing plane of the LISA detector. 

Figure 6. The same as Fig. 5, but for the non-preeessing OMEGA detector. The entire 

gravitational wave power is now concentrated in a single peak. 

Figure 7. The all-sky map of peak spectral power in the demodulated signal for one-degree 
pixels. The signal being analyzed represents the received LISA time series for three binary 
gravitational wave sources plus random noise at the predicted LISA level. The parameters of 
the three sources are: SOURCE 1 — /i = 3.3 x lO'^^, / = 3.0 x lO'^ Hz, 9 = 60° , (i) = 120°, 
i = 10°, i/j = 84°, (Po = 0; SOURCE 2 — /i = 4.0 x 10-^2, / = 4.0 x lO'^Hz, = 78°, 
= 20°, i = 40°, V = 0°, ^/)o = 0; and SOURCE 3 — h = 3.7 x lO'^^, / = 3.5 x lO'^ Hz, 
= 12°, (f) = 280°, i = 0°,il; = 84°, (f)o = 0. 

Figure 8. The same type of all-sky map for the same signals as Fig. 7, but for the non- 
preeessing OMEGA detector. The relative weakness of the power near the ecliptic (the 
= 78°, (j) = 20° location) is due to the non-uniform sensitivity of the OMEGA detector - 
more sensitive than LISA near the ecliptic pole and less sensitive near the ecliptic plane. 
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